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A sea, swell and surf program is improved, tested and evaluated on a 
micro-computer (HP-9845B). Sea swell is calculated by a two dimensional 
spectral model. The energy balance equation is tested for different cases of 
wind velocities and water depths. Satisfactory agreement is observed 
between the offshore model and expected wave heights for a 15 knot wind, 
but the model overbuilds wave energy for a 30 knot wind. Wave 

transformation is described by a one dimensional random wave model in which 
the wave heights are described using the Rayleigh distribution. The obtained 
solution of the random wave field is used to predict the longshore currents. 
An empirical formula for determining the breaker parameters is developed, 
based on beach slope and incident wave steepness. The improved model is 
tested using an undulated bathymetry to validate the model physics. The 
model outputs of wave height and current are compared with data acquired 
from a wave tank and natural beaches. The model is found to accurately 
forecast wave heights, breaker location, breaker type and longshore currents 
for several sets of conditions. Model limitations are discussed and 
recommendations for further improvement are made. 
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I. INTRODUCTION 



A. OPERATIONAL REQUIREMENT FOR WAVE FORECASTING 

Wave and surf forecasting has been an important part of Naval 
Amphibious operations since the beach landings of World War II. The need 
for the prediction of ocean and surf waves for military purposes led to the 
first real attempt to quantitatively model ocean and beach waves. Ocean 
engineers and physical oceanographers, such as M.P. O’Brien, H.U. Sverdrup 
and W.H. Munk, made great strides in sea, swell and surf modeling during 
the 1940 - 1955 time period (Bascom, 1980). The wave and surf predictions 
made by the U.S. Navy are still based on these works. During the 
intervening time period, significant advances in knowledge of ocean wave 
dynamics have occurred, and it is time to update prediction techniques. 

Although warfare technology has advahced along all fronts, the standard 
amphibious beach assault is stUl a viable and preferable method of projecting 
forces ashore for many tactical scenarios. Hqwever, amphibious assault craft 
are vulnerable to many types of wave dynamics. In general, landings are 
affected by the initial ocean sea and swell where the boats are launched as 
well as by the breaking waves in the surf zone, along with the accompanying 
longshore currents. Since in a major landing tens of thousands of lives are 
at risk, an accurate assessment of the sea, swell and surf is vital to a 
successful amphibious landing operation (Joint Surf Manual, 1976). 

In addition to the assault landing itself, an amphibious operation normally 
concludes with the construction of a temporary harbor facility to 
accommodate the resupply of the troops ashore. The survivability and design 
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characteristics of these structures depend in large part on the environment 
around them. The ability to model wave height and direction, and longshore 
currents is necessary to operational planners and coastal engineers in 
designing and emplacing these facilities • 

% 

In the past twenty years, separate models have been developed to deal 
with ocean waves, breaking surf and longshore currents. In 1981, the U.S. 
Navy contracted to have a model developed which would combine three 
complementary models. The goal of this comprehensive Sea, Swell and Surf 
Program (SSSP) model was to predict the wind generated wave height, 
calculate the surf zone characteristics based on the wind generated energy 
propagation, and use a two dimensional grid to compute the longshore current 
field. The model, developed by Wang and Chen (1983), was written in 
FORTRAN and designed to run on a main frame computer. To facilitate the 
model’s use by Naval Oceanographers in the fleet, the model was converted 
into BASIC and implemented on the HP-9845B/275 micro-computer 

(Devendorf, 1985). The primary goal of this thesis is to modify the SSSP 

model to improve its computational efficiency. 

B. WAVE THEORY BACKGROUND 

The water surface elevation in the ocean can be viewed as a large 
number of individual wavelets at various frequencies. The waves that con- 
tribute the most energy to the ocean wave system are wind generated waves 
with periods of 1 to 30 seconds. The size and period of these waves are a 
function of the velocity of the wind, the duration of the storm and the 
distance over which the wind blows (fetch). 
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The simplest theory is linear, small amplitude wave theory. Sea water 
is assumed to be incompressible and homogeneous, and surface tension forces 
are assumed to be negligible , The solutions of this theory describe the 
behavior of the individual components of ocean wind waves (Phillips, 1977), 
As waves enter into shallow water, wave speed varies with the local water 
depth and refraction occurs. 

Initially, refraction of monochromatic waves was modeled using wave ray 
theory of a single wave (Arthur, 1949), It was not until the mid 1960’s 
that the development of fast computers and the implementation of better 
physics enabled numerical analysts to study the generation of wind waves 
over an entire frequency spectrum. Early researchers in the wave spectral 
analysis field introduced direction variables to study straight and parallel 
contour wave spectrum transformation (Karlsson, 1969 ) and included the 
effects of bottom dissipation and local wind generation. 

As computers advanced, more complex refraction models were developed, 
Noda et al. (1974), using a relaxation finite differencing scheme, solved for 
the stationary wave spectral transformation. The model that is under study 
by this thesis was originally developed by Wang and Yang (1981) and 
included the effects of bottom friction, Chen and Wang (1983) improved the 
model by including the effects of non-stationary waves, 

C, CURRENT MODELING TECHNIQUES 

Two numerical modeling techniques are currently used in wave and surf 
forecasting applications. The finite element method, while highly complex, 
has the advantage of being able to incorporate irregular localized bathymetry 
and has stable non-linear convergence properties, Wu and Liu (1985) 
introduced the method in their study of wave induced nearshore currents. 
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The finite difference method is a widely used numerical method. It is 



based on a linear differencing approach. Finite differencing normally uses a 
rectangular computational grid. This method was used by Noda et al. (1974 ) 
and Shaiu and Wang (1977) in their studies of nearshore circulation and wave 
energy transformation, respectively. The method is easier to implement than 
the finite element method due to its use of a regular grid describing the 
bathymetry, and is employed in the development of the SSSP (Wang and 
Chen, 1983). 

D. OBJECTIVES 

The objective of this thesis is to improve a sea, swell and surf program 
which is implemented on the HP 9845-B/275 micro-computer. The program is 
divided into subroutines, which include the open water wave generation 
module , wave transformation , the surf zone breaker calculations and the 
longshore current routines. To test the model, the output from the improved 
model is compared with data bases acquired from two natural beaches as well 
as several wave tank experiments. A users manual has been prepared and an 
attempt has been made to interface graphic output with the current model. 
Limitations of the model are discussed and future improvements suggested. 
This model is proposed for use by operational planners and Naval Ocean- 
ographers at sea . 
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II. OVERVIEW OF THE SSSP MODEL 



A. OVERALL PROGRAM DESIGN 

The Sea, Swell and Surf Program (SSSP) was developed by Wang and 
Chen (1983). The SSSP is essentially three complementary models which are 
merged to calculate offshore wave height fields, surf zone information and 
longshore current velocities. The three modules, ocean, surf and longshore 
currents, can be run independently or the offshore wave energy can be used 
as input to the surf and current models. The SSSP is a numerical model 
which uses finite differencing to solve the governing equations. 

The original version, as written by Wang and Chen, was written in 
FORTRAN and designed to run on a mainframe computer. The version 
discussed in this thesis was converted to BASIC and implemented on an 16 
bit HP 9845-B/275 micro-computer (Devendorf, 1985 ). The micro-computer 
version consists of a control program with subroutines to calculate wave 
direction, height, and number, and output surf, swell and current 

information . 

The SSSP considers an offshore wind generated energy spectrum which is 
used to estimate the wave growth based on initial inputs of wind speed and 
wave height or spectral energy bands. As the waves traverse down the user 
defined grid, changes in the deep water ‘ energy spectrum due to bottom 
friction and irregular bottom changes, are taken into account. If the surf 
zone model is to be run, a significant wave height , calculated from the area 
under the wave spectrum, and the peak frequency are used to provide a 
monochromatic wave height input into the surf model. It is also possible to 
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input monochromatic wave data directly into the surf module to compute the 
surf zone characteristics without running the computationally intensive 
offshore module (Wang and Chen, 1983). 



B. SSSP FLOWCHART 

The SSSP is controlled by one main program which, in turn, can call any 
of twelve subroutines. A general flowchart of the entire SSSP shows that 
the SSSP is divided into two main modules (Fig. 2.1). The user is prompted 
for all required information in a menu-driven, interactive fashion. Initially, 
the user chooses whether to run the offshore or the surf model. General and 
specific input is requested for the module being run. The open water or surf 
condition modules are calculated and the proper output is generated by either 
the OUTSW or OUTSRF routines. 

A more detailed flowchart of the open water condition module is 
presented as Figure 2.2. OPEN computes the wavenumber as a function of 
depth and then loops through a calculation series for each energy band. An 
adjustment is made for the angle of the swell/ waves with respect to the grid 
axis. The DIRECT and HEIGH2 subroutines are called and the energy from 
each spectral band is summed and converted to wave height. 

The original SURFCON subroutine (Figure 2.3) used the most energetic 
energy band computed by the OPEN module as its frequency input, or used 
data entered directly by the user. The wavenumber was computed, the 
DIRECT and HEIGH2 subroutines were called, the average wave height was 
calculated and a series of breaker height and locations were predicted. 

If longshore currents were desired, the NEAR Cl R subroutine was invoked. 
The model, as implemented on the HP-9845 , used an iterative process on a 
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General flowcharT of The Sea, Swell and Surf Program 
(SSSP). 
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Figure 2.2. Flowcharts of the Ooen water, wave Direction anfl wave 
Heiaht subroutines (OPEN. DIRECT and HE1GH2). 
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two dimensional grid to solve for the mean horizontal velocities. It was 
found that this scheme was computationally unstable as well as requiring as 
long as twenty-three hours to converge to a solution (Devendorf, 1985 ). 

To take advantage of the most recent advances in breaker height 

modeling, and increase the computational speed, the SURF CON and its 

ancillary subroutines were revised to use a one dimensional random wave 
model. detailed flowchart of the new SURFCON module (Fig. 2.4), which 
replaces the original, shows the simplifications implemented by the new 
model. DIRECT2 is a simple SnelUs law calculation, applicable over straight 
and parallel contours. HEIGH3 uses a random wave model which incorporates 
the probability density function for breaking waves. The output of the new 
SURFCON module includes breaker height, number of breakers, wave direction 
and surf zone width. 

LONSHOR is a new longshore current subroutine (Figure 2.4), which was 
implemented in place of the unstable NEARCIR to simplify and speed the 
calculation of the longshore currents. This model uses a simple radiation 

stress balance to calculate the current parallel to the beach. The next 

section will discuss the physics of the SSSP modules in detail. 

C. PROGRAM MODULES 

1. Offshore Module 

When the SSSP is run, a menu of initial choices is presented. If the 
user chooses the "OPEN WATER" option, the OPEN subroutine is invoked. 
This offshore module computes the wave height and direction fields as a 
function of water depth, wind speed and direction, fetch, and the initial 
energy field. The initial energy can be described by a single significant 

wave height, direction and period, or as an implicitly entered energy 
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Figure 2.4. Flowcharts of NEWSURFCON, D1RECT2, HE1GH3 and LONSHOR 
subroutines used in the new model. 
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spectrum . 



Water depth for the matrix can be read in from a data file or 



entered directly from the keyboard. A tidal height correction constant can 
be added to the depth field to take into account the state of the tide . 
The bottom type can be specified, which determines the coefficient of friction 
for the wave height calculations, 
a. Wave Refraction 

Once all the required input is entered, the program begins its 
calculation cycle. The first field that is calculated is the wave direction 
matrix, using a steady state conservation of wavenumber (Phillips, 1977 ): 

cos6) - T^(k sin6) = 0 (1) 

o A d Y 

where : 

. , 2tt . “1, 

k = wavenumber = (m ) 

±j 

L =• wavelength (m) 

6 = wave direction (clockwise angle from 

the X axis ro the wave ray (see 
Figure 7 ) ) 

X =. offshore grid spacing (user specified 
(see Figure 5) ) 

y = longshore grid spacing (user specified 
(see Figure 5) ) 



The wavenumber is calculated, for the frequency band being considered, using 
the linear dispersion relation: 

0 = (g k tanh (k h))^^" (2) 
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where : 



c = angular frequency = 2zf (sec ) 
g = acceleration of gravity (m sec ) 
h = water depth (rn) 
f = frequency = ^ (sec 

T = wave period (sec) 

Since (2) is transcendental in k, a sixth order polynomial fit of (2) is used 
(Hunt, 1979 ) to provide a faster method of computing k than the iterative 
Newton’s method used in the original Wang and Chen (1983 ) version of the 
SSSP. The wavenumber k is given as: 

k = ^-r-jylS + (1 + 0.66667S -f 0. 355503^ 

(gh)*^^^ 

^ 0. 160845^ + 0.06320s'^ ^ 0.021743^ 

+ 0.006543^ + 0.001713"^ + 0.000393® 

9 -1 1/2 

+ 0.000113 ) (3) 



where : 



3 






Once the wavenumber is determined, subroutine DIRECT is 
called to compute the wave refraction of the wave components over the 
specified grid. DIRECT starts with initial values of the wave angles for each 
grid point, and uses a finite differencing scheme developed by Noda (1972) 
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to solve the differential equation for . The differential equation is 

center differenced in the offshore ( X ) direction and forward differenced in 
the longshore (Y) direction* The finite difference scheme has greater weight 
on the forward differencing term to increase computational stability (Abbott, 
1979 ) and is given as: 



6 . . = s in ^ { r— ^ — [T(ksin0. , ^_,)-(l-2T)(ksin9._^, ) 

^ f J 



+ T(k sine.^^^ 



- ('k cos 6 • • 1 ) 3 } 

1 , J - ± 



(4) 



where : 



T = weighting factor = 0.25 

A first guess for the wave angle is calculated using Snell's law 
solution for straight and parallel bottom contours. The numeric scheme is 
iterated until £ . ^ is within 0.005 of the previous calculation. The 10 X 13 
grid point system used by the SSSP is described by Figure 2.5. The initial 
differencing scheme used by the model is a simple forw'ard stepping scheme 
(Fig. 2.6). The convention for the angle is illustrated in Figure 2.7. 
Theta is measured with respect to a ray drawn perpendicular to the beach. 
Therefore, a normally incident wave would have a Theta value of 180^. 
After the direction matrix has been specified, the subroutine HEIGH2 is 
called to calculate the height transformation as the wave moves down the 
grid (towards the shoreline). 
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Shoreline , 
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Figure 2. 6, Numerical differencing scheme used by the OPEN water 
module. 
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Figure 2.7. Theta convention used by the model. 
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2 . Wave Height Transformation 



This subroutine takes into account the wave direction (previously 
calculated), the wind generation, and the bottom dissipation, HEIGH2 
employs the Phillips-Miles growth mechanism for the wind generation 
calculations. An initial energy field is developed for each grid point by 
assuming that energy is conserved over straight and parallel contours: 

■^(E Cc cos 6 ) = 0 (5) 

where : 

E = wave energy (joules) 

Cg = group velocity (m sec*^) 



The initial energy field is derived by integrating (5) to yield: 

Cg cos 0 

E = E — - ^ 

o Cg cos 6 



where the subscript (o) is the initial condition at the offshore boundary. 

To simplify and increase the model’s efficiency, Chen and Wang (1983 ) make 
four assumptions: 

a. Winds are steady 

b. Wave energy within a spectral band is restricted to that 
band (Allows superposition of wavelets). 



24 



c. Each frequency component can be described by a single 
mean direction, G . 

d. No wave-current interactions. 

These four assumptions allow the steady state energy flux for each spectral 
energy band to be written as ( Longuet-Higgins and Stewart, 1960; 1962 ): 



jg 



— ( S ( f ) Cg cos - ) 

OA 



>g — (S (f ) Cg sin t) 
-^0 1 



"d 



f) * 



(f ) 



ill 



where : 



S (f ) 



0 




spectral energy flux 

de.nsity of fluid ()<gm '') 

energy dissipation due to 
botrom friction 

energy generation due to 
wind stress 

mean wave direction 



The wind generation term £ ^ is the sum of the linear and exponential 
contributions : 



£ 

S 




( 8 ) 
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wnere ; 



P 



Phillips generation 



pg:x 



= Miles generation = pg$[S(f)] 



The energy generation term is based on the Phillips-Miies growth mechanisms. 
Starting with a flat ocean, wave growth is initially linear until larger waves 
are present, after which the growth is exponential. The linear growth term 
(Phillips, 1957) is based on random atmospheric pressure variations acting on 
a perturbed sea surface. The SSSP first calculates the pressure fluctuation 
term as a function of wavenumber and frequency: 



P(k, c ) 



where : 



- ^ 6 

6.13 '10 "K 



^ 7 



! k sin c ) 




'' 1 



(k cos 




(9) 



W = 'wind soeed in knots 



W 



1 — 1 

> 


= 


0 .3 3/1 


. 28 


“I 


= 


0 . 52/1 


.95 






angle 


between wind and 



After P(k, 0 ) is calculated, the linear growth term, a , is determined by: 



a 






p g 



P ( k , G ) 



( 10 ) 



26 



where : 



0 ^ = density of water (kg m 

The exponential growth term becomes important after waves have formed 
(LaBlond and iMysak, 1978). When flow separation occurs in the lee of a 
wave crest , momentum is transferred to the wave because of the pressure 
differential resulting in exponential wave growth (Miles, 1957; 1959a, b; 

1962) : 

6 = (^) [ (p)COS 6 -b] (11; 

Z h L- 

where : 

a = constant = 5.0 

S = ratio of air to water density 

C = wave celerity (m sec 

b = constant = 0.90 

The Phillips-Miles wind wave generation has no mechanism for 
shutting off the wave grow'th. In nature, wind waves will grow as a 
function of wind speed, duration and fetch. But after growing to a certain 
size and steepness, the waves will become unstable and break. Therefore, it 
is necessary to have a limiting function so that the waves will only grow to 
a ’’fully arisen" sea. Several functions have been proposed to introduce an 
"equilibrium value", above which the wave growth stops. The SSSP uses a 
Pierson and Moskowitz (1964) fully arisen sea in deep water. The high 
frequency tail is described by a Phillips (1957 ) equilibrium condition. 
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modified by Kitaigorodskii et al. (1975 ) and Thornton ( 1977 ) for intermediate 



or shallow water; 



where : 



S(f,6) = B ( c) r (V7j^) 



( 12 ) 



H(6) 



r(W^) 



w 



n 



8 4 

= -j 71 cos 6 accounts for angular 
spreading 



,-2 



<'"n> 11 



2 w-i : (w ) 

n n 

sinh[2W^Z(W ] 
n n 



-1 



9 



Z(W ) meets the condition that I tanh( W '^ ) = 1 
n n 

r (W ) = 1 in deep water and for W < 1, T(W ) — 1/2 W ^ and S(f, c ) = 
n ^ n ’ n n ’ 

“3 “2 

1/2 q g h df H( C ) where q is on the order of 10 (SSSP uses a value of 

0.073) . 



Bottom dissipation, is calculated as the work done on the 

bottom due to the bed shear stress. Assuming the usual quadratic bed shear 
stress law, the wave energy dissipation for a particular spectral frequency 
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band is described by Hasselmann and Collins (1968 ): 



where : 



c^(£) 



!“bl 



( 13 ) 



'f 

U. 



= bed shear stress coefficient = 0.01 



= total flow field bottom velocity 



Ufalf) 



wave induced water particle bottom 
velocity for a specific spectral 
component 



Using linear wave theory transfer functions to relate the surface elevation 
energy spectrum to the bottom velocity field, the frequency dependent 
dissipation function £^(f) can be expressed as: 

C^(f) = P CjT(f) S(f) (14) 



2 2 2 2 

where the transfer function T(f) = [ g k / o cosh (kh)). The total wave 
induced bottom velocity 1 I is calculated by integrating the spectral 
transfer function across the entire frequency spectrum: 

1/2 

= [ / T(f) S(f) df] (15) 

0 

Bottom dissipation in the SSSP is only calculated when deaLng with 
intermediate or shallow water waves. If deep water wave heights are being 
computed, the value ^ ^ is set equal to zero (Wang and Chen, 1983). 



29 



3 • Longshore Currents 



One of the most serious problems with the original model was the 
instability of the nearshore circulation calculations. As much as 23 hours of 
CPU time were required for the subroutine to converge to an answer. It 
was decided that a one dimensional longshore current model which assumes 
straight and parallel bottom contours would be substituted for the unstable 
two-dimensional model. This model uses a simple alongshore radiation stress 
balance to calculate the mean longshore current. Although no circulation 
information is derived, the new longshore current subroutine, LONSHOR, 
provides the longshore current velocity efficiently and accurately. The 
physics and details of the LONSHOR subroutine are outlined in Chapter II. 

4 . Surf Zone 

The original SURF CON module made calls to the same wave direction 
(DIRECT) and wave height (HEIGH2) subroutines as the offshore module. 
One disadvantage of this scheme is that these subroutines are very 
computationally intensive, especially when considering the narrow width 
between grid points. A faster wave height model was needed. 

The state of the art in breaker height models uses a probability 
density function to describe breaking wave heights (Thornton and Guza, 
1983). Some limitations are imposed by using this new model, however. 
The model assumes that the beach has a local straight and parallel bottom 
contour. This bathymetric restriction allows longshore features such as 
sandbars to be modeled accurately, but offshore feature such as channels or 
canyons are smoothed out by a depth averaging subroutine. It is felt that 
this assumption is reasonable , in light of the fact that many beaches of 
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operational interest are of the "straight and parallel" type within the first 
two hundred meters from shore. 

The above assumption means that a simple Snell*s law calculation of 
the wave direction is possible. A short subroutine, called DIRECT 2, is 
implemented to calculate the wave direction matrix. The use of the new 
height and direction modules, HEIGH3 and DIRECT2 have increased the 
efficiency of the SURFCON calculations immensely. The average run time for 
the surf zone model predicted wave heights has decreased from over four 
minutes to less than 30 seconds. 



31 



III. MODEL IMPROVEMENTS 



A. OFFSHORE WAVE HEIGHT GROWTH 

The current version of the SSSP, converted to run in HP BASIC on the 
HP-9845B/275 micro-computer, gives reasonable values for wave height 
growth for wind speeds of 15 knots or less (Devendorf, 1985). However, for 
higher wind speeds of 30 knots, the SSSP exhibits a growth rate that is 
significantly faster than the growth of the JONSWAP spectrum (Bishop, 
1983), used for comparison in this analysis. The JONSWAP spectrum was 
chosen as a basis for comparison because the JONSWAP curves were 
experimentally derived from the limited fetch North Sea area. It is foreseen 
that the SSSP model may be applied to similar operational areas. 

For 30 knot winds, the model^s wave growth not only occurs too 
quickly, reaching a maximum value within 50 nautical miles, but the 
magnitude of the growth is too high (maximum wave height of 6.0 meters 
compared to a JONSWAP maximum of 4.2 meters. Before sensitivity tests 
were run on the wind generation portion of the model, the code was checked 
carefully to ensure that the equations were being implemented correctly. 

To determine the sensitivity of the model wave growth, several runs 
were made with changes to the key growth parameters. These parameters 
include A, the Phillips linear growth term. Equation 9 (Phillips, 1957); B, 
the Miles exponential growth term. Equation 10 (Miles, 1957; 1959a, b; 

1962), and R, a weighting factor set equal to 1.0 in the original model. 
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The first run was used as a standard to compare the modifications 
against. The inputs for the OPEN module were: 

Grid Spacing X = Y = 10 NM 

Wave Direction Theta = 180° (normal incidence) 

Wind (W) = 30 and 15 kts @ 270° 

Shoreline = 270^ 

Bandwidth (D£) = 0.01 Hz 

Frequency bands = 0.05 - 0.14 Hz 

The output for these "standard” model runs are shown as Figures 3.1 and 
3.2. The overbuilding of the open water waves for W = 30 kts is evident 
when compared with the JONSWAP spectrum wave growth (Fig. 3.1). For a 

more moderate case where W = 15 kts, the model wave growth shows 
satisfactory agreement with the JONSWAP growth (Figure 3.2). To see the 
impact of changing the weighting factor, R, two runs were made setting R 
equal to 0.5 and 0.1 respectively. The predicted wave height with R equal 
to 1, 0.5 and 0.1. is shown in Figure 3.3. Setting R equal to 0.1 (chain 
dashed line) gives a good fit for high wind speeds but underbuilds the wave 
field for lighter winds (Fig. 3.4). Clearly, a simple change of the 

weighting factor R will not resolve the problem. 

The Miles exponential growth term, B, is in the denominator of the 
energy calculation so its value was increased to determine the growth 
retarding effect it would have. Two runs were made with B doubled (chain 
dotted line) and with B increased by a factor of 10 (not shown) , depicted 
by Figure 3.5 The result of increasing the exponential growth term is to 
flatten out the upper portion of the growth curve which is not the correction 
that needs to be applied. 
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FE'l'CIl DISTANCE (NM) 

Figure 3,1. Model output (solid) versus JONSWAP curve (dastied) for W 
; = 30 kts. Tfie SSSP substantially over develops wind 
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Figure 3.2. Model output (solid) versus JONSWAP curve (dnsfied) for 
= 13 kts. For low wind speeds, there is good agreement 
between JONSWAP and SSSP. 
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l-ETCH DISTANCE (NM) 

Figure 3.3. Sensitivity test; R = 0.5 (chain-dot) and 0.1 (chain-dash) 
for W = 30 kts. R is a weighting term which is set equal 
to 1,0 in the original SSSP model. The unclianged SSSP 
model (solid) and JONSWAP (dashed) are also shown. 
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Figure 2 . 4 . Sensitivity test; R = O.l (chain-dot) for W = 15 kts. 

Since each is rnul tipiieci by R, a reduced R vaiiie will decrease 
the low wind speed case as well as the iugh, Ttie uncimnged 
SSSP modei (solid) and JONSWAI* (daslied) are also shown. 



Because the model’s rapid growth appears to occur in the linear region 
of the growth curve, the Phillips linear term, A, was analyzed closely to 
see if a simple constant multiple would bring the wave growth in line with 
the JONS WAP curves. A constant multiple of 0.125 is applied to the A term 
(Fig. 3.6). This method gives a better fit to the JONSWAP curve for high 
wind speeds, although the SSSP still over builds the wave heights. To test 
this multiple at lower wind speeds, a model run was made for a wind speed 
of 15 knots, for a bandwidth of 0.05 Hz and frequencies of 0.05, 0.10, 

0.15, and 0.20. The results of the model run are compared with the 
equivalent JONSWAP curve and the curve from the existing model (Fig. 3.7). 
The fit, while not exact, is reasonable, suggesting that the A term is 
responsible for the overbuilding of the model wave heights. 

The actual cause of the overbuilding of wind waves by the SSSP is 
unknown. Several ad hoc "fixes” to the model are suggested with the 
recommendation that further research on this problem be conducted. A 
simple constant multiple of the linear growth term, on the order of 0.10 to 
0.15 is one recommended empirically derived modification to the wind 
generation calculation. Another possible modification is to have the 
weighting function, R, decrease as the wind increases, to act as a damper 
on the over generation. At present, the SSSP wave generation routine 

remains in its original form. It is felt that the model is adequate for low 
wind speed generation. Operational planners should be cautioned about the 
high wind speed generation problems. 

B. WAVE REFRACTION MODULE 

To test the SSSP model’s stability during the wave height prediction in 
the coastal zone, a sinusoidally varying bottom wais explicitly entered using 
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Figure 3.5. Sensitivity test; R = 2R for W = 30 kts ( chnin-dot ) . B 
is the exponential growth term. Increasing the R value t>nly 
tends to flatten out the upper part of the curve. The unchanged 
SSSP model (solid) and JONSWAP (dashed) are also shown. 
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FETCH DISTANCE (NM) 

Figure 3.6. Sensitivity test; A term reduced by constant multiple of 
0.125 (chain-dot) for W = 30 kts. The unchanged SSSP 
model (solid) and JONSWAP (dashed) are also shown. 
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Figure 3.7. Sensitivity test; A = O.I2f) A ( chain-dot ) for W = lf> kts. 

The unclianged SSSP model (solid) and JONSWAP (dashed) 
are also shown. 



the following equation: 



h, , = 0.025X + ch, {^){1 - (r^))cos(^) (16) 

(x,y) b 

where : 



E = bottom perturbation term (0.2) 

h, = breaker depth (2.5 m) 

^ b 

X, = offshore distance to breaker line 
(100 m) 

= wave length of bottom perturbation (240 m) 

This equation provides a smoothly varying bottom profile in both the X 
(offshore) and Y (longshore) directions. To check that the bottom profile 
was being accepted correctly by the model, and that the calculated results 
were reasonable, approximately 25 model runs were completed with varying 
input parameters. Selected cases are presented herein for discussion. In 
the cases presented, the offshore waves are normally incident to the 
shoreline. Obliquely incident model runs were made with Theta ranging from 
+10 to -10 degrees from normal, for comparison with the normally incident 
case . 

1 . Model Experiments 

The input parameters in the model are: significant wave height, 

(1 m); wave period, T (12 s) ; offshore boundary depth, h(l2 m); artificial 
viscosity, T (0.0); and bottom perturbation term, £: (0.2). The beach 
faces 270° and the normally incident wave rays come from 270°. This gives 
a Theta value of 180° according to the model convention shown in Figure 
2.7. 
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The model results for wave heights and directions are shown in 
Figures 3.8-3.15. The longshore depth profiles for locations inside (X=3) 
and outside (X=8) the surf zone are chosen to examine the wave refraction. 
Grid point X = 3 implies a distance of 40 meters offshore, well inside the 
breaker line, while X = 8 is 14 0 meters offshore and outside the breaker 
line. The depth profile presented in Figure 3.8 shows that the depth varies 
sinusoidally in the longshore direction . The beach profile along grid point Y 
= 3 is shown as Figure 3.9. As shown, the depths are slightly deeper than 
the plane beach inside 100 meters (breaker position), and shallower outside 
100 meters. 

As waves enter the surf zone, their increase in height is a function 
of water depth. For the points outside the breakers, the largest waves are 
expected to occur where the water depth is the shallowest due to shoaling. 
As can be seen from Figure 3.10, the peak waves (dashed) are towards the 
edges of the domain, but they do not coincide with the minimum in depths at 
0 and 240 meters. This anomaly is a result of the boundary conditions of 

the model. Also note the asymmetric peaks at 40 and 200 meters. This 

implies errors in wave height prediction due, in part, to the wave direction 
calculations. These errors will be addressed later in this discussion. 

Wave direction is calculated in terms of Theta. An angle greater 

than 180 tends to turn the wave ray to the left while Theta less than 18 0 

turns the ray to the right (See arrows in Figure 3.11). The model at least 
quantitatively predicts the wave refraction pattern correctly for all cases, 
consistently turning the rays to "attack” the shallower water depths outside 
the breaker line. The models boundary conditions force ^ ” ^2 ^ n-1 

= 9^ where G is the wave direction at the n^^ long shore grid point. This 
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LONGSHORE DEPTH PROFILE (EPS = 0.02) 
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LONGSHORE DISTANCE (M) 

Figure 3.8. Longshore depth profile for X = 3 (s{)lid) and X 
: (dashed). The sinusoidally varying bottom profile 

used as a depth field for testing the SSSP. 
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Figure 3.9. Offshore (feplli profile for Y = 3 (solid). Depths me 
(lee|)er/ shallower than a plane beach profile of slope 0.025 
(dashed) . 
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Figure 3.10. Wave Heights at X = 3 (solid) and X = 8 (dashed). Note 
asymmetric peaks at 40 and 200 meters. Wave heights 
inside the breaker line are related to the bottom depth 
only . 



WAVE DIRECTION FIELD (THETA = 100) 
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Wave DirocLions at X = 3 (soli(i) and X = 8 (dashed). 
Arrows indicaLe the direcLions the wave rays point witli 
respect to a ray normal to tlie beach. Dotted tines at 



accounts for the departures of the wave direction at the edges of the 
domain. The dotted lines are the interpretation of how the wave direction 
patterns should behave. In the middle of the domain, where both cases 
exhibit a nodal point, the wave direction passes through 180 or normal to 
the beach as expected. 

The final test conducted compared the effect of Tau on the wave 
direction calculations. Tau is an artificial viscosity parameter which is 
intended to be a smoothing term for the calculation of wave angle. 
Changing Tau from 0.00 to 0.25 had a negligible effect on the wave height 
calculations, but the difference in direction increased with decreasing 
distance from shore. The case where X = 8 is plotted in Figure 3.12. The 
introduction of the Tau term (dashed) decreases the magnitude of the wave 
angle, but the maximum (minimum) values are not exactly in phase with the 
Tau = 0.00 case. The maximums are shifted slightly toward the edges of the 
domain. This is the expected behavior of the smoothing term. 

The dissipative model exhibits unexplained behavior inside the 
breaker line as depicted by Figures 3.13-3.15. AtX = 5 (Fig. 3.13), 
the Tau model reduces the magnitude of the wave direction field. At X = 4 
(Fig. 3.14), spurious asymmetry appears in the direction calculations with 
the Tau term applied. Furthermore, at X = 3 (Fig. 3.15), the effect of Tau 
is to increase the magnitude of the wave direction. This apparent angular 
deviation has no effect on the model’s wave height calculations because , 
inside the breaker line, wave height becomes a function of depth only. 
However, this instability does effect the longshore current calculations, 
which are sensitive to local wave direction inside the surf zone. 
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Figure 3.12. Wave F)lrections at X = 8 for Tail = 0.0 (solid) and 0.2ri 
(dashed). The condition tliat 0^ = (j ^ is a function of 
tlie original boundary conditions. 
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LONGSHORE DISTANCE (M) 

Figure 3.13. Wave Directions at X = 5 for Tan = 0.0 (solid) and 0.25 
(dashed). The expected modifying behavior of Tan is 
exliibited . 



WAVE DIRECTION FIELD (THETA = 100) 
TAU = 0.0 AND 0.25 AT X = 4 
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LONGSHORE DISTANCE (M) 

Figure 3.1 A. Wave DirecLions aL X = 4 for Tau = 0.0 (solid) and 0.25 
(dashed). As the wave approaclies the shoreline, the 
'I’au rnodificaLion shows irregular behavior. 



WAVE DIRECTION FIELD (THETA = 180) 

TAU = 0.0 AND 0.25 AT X = 3 
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LONGSHORE DISTANCE (M) 

Figure 3.15. Wave Directions at X = 3 for Tau = 0.0 (solid) and 0.25 
(dashed). Contrary to its intent, close to the shoreline, 
Tau actually increases the magnitude of the wave 



Another indication of a numeric error in the wave direction 
calculation was noted in the exact center (Y = 7) of the domain. Since the 
depth profiles exhibit a maximum (minimum) value at the domain center, 
theory suggests that the waves remain normally incident at these nodal 
points. Therefore, the Theta value should be 180 degrees at the center of 
the domain. The model output, both with and without Tau, exhibits a 

V 

monotonically increasing error as seen in Figure 3.16. 

2. Improvement of the Numerical Scheme 

In an effort to eliminate the numerical instability in the wave 
direction module , the finite differencing scheme was analyzed to see if an 
improvement could be made. The original scheme, pictured in Figure 2.6, 
uses a simple centered suid forward differencing method and is fairly robust in 
the longshore (Y) direction but is less so in the offshore (X) direction. To 
increase the accuracy of the numerical scheme, the method pictured in 
Figure 3.17 was substituted in the DIRECT subroutine for the original 
scheme. This particular method is a well accepted technique which has 
proven to be stable and free of damping (Haltiner and Williams, 1980). 

Several model runs were made to test the new technique. Using the 
standard input parameters euid setting Tau equal to 0, the first result was a 
gain in accuracy of the wave direction calculations. The domain’s central 
directions were exactly 180^ (to three significant figures), as seen in Figure 
3.18, and the refraction angles were symmetric about the center of the 

domain as expected due to the symmetrically varying bathymetry (Figure 

3.19). 

A second result of the implementation of the improved numerical 

technique is the stabilizing effect on the wave height calculations. The 
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WAVE DIRECTION FIELD (THETA = 100) 

MONOTONIC INCREASE OF THETA 
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Figure 3.16. Monotonic increase of Theta (solid) . versus the expected 
value of 180° at beach center. 
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A more robust differencing scheme was subsLitutecl for the 
original numeric scheme seen in Figure 2.6. 
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OFFSHORE DIS'i’ANCE (M) 

Figure 3.18. Increased stability in Theta calculations resulting from the 
new numeric scheme (com 0 are witli Figure 24). The 
expected value of 180.000 is achieved with the new 
scheme . 
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Wave Directions for Tan = 0.0 at X = 8 (soliil) usin^ tlie 
now numeric sclieine compared willi the old scheme 
(dashed). Doited lines show ex[)ected behavior at the 
boundaries. 



previous asymmetric wave heights are now symmetric with respect to height 
and peak location (Fig, 3.20). 

3. Boundary Conditions 

The SSSP model imposes no flux boundary conditions on the 
refraction and wave height calculations so that the left and right boundaries 
are specified as follows; 

Q(M,1) = Q(M,2) 

Q(M,N) = Q(M,N-1) 

where: Q = the quantity being calculated 

M = offshore row number (10 to l) 

N = the last longshore column (13) 

The offshore row is completely specified by the initial conditions. The 
calculations then proceed from offshore to onshore and from left to right on 
the grid. 

No flux boundary conditions account for the unusual behavior of the 
wave height and direction curves at the edges of the domain. To satisfy the 
condition of well-posedness for the applied model, cyclic boundary conditions 
are substituted such that: 

Q(M,1) = Q(M,N“1) (17) 

Q(M,N) = Q(M,2) 

A plot of wave height versus longshore distance at X = 8 using the same 
inputs as the previous case is presented in Figure 3.21. Note the shift of 
the peaks toward the edges of the domain as well as the abrupt increase in 
the magnitude of the wave height. These numerical results now match quite 
well with theory which requires that the wave rays turn toward the 
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shallowest water, implying that the largest wave heights should occur at the 
edges of the domain. 



The wave direction calculations also show improvement. Close 
agreement between the cyclic boundary condition values and the expected 
values (dotted) is shown in Figure 3.22. Additionally, the average value on 
the boundary segment agrees closely with the exact value (solid circles). 

4 . Results 

Based on the results of many model experiments, the new model 
appears to offer an increase in stability and computational accuracy. The 
improved wave direction subroutine consistently turns the wave rays in the 
correct direction . The magnitude of the direction changes is qualitatively 
correct. The wave height subroutine works well both inside and outside the 
surf zone. Wave heights are highest at the shoal and lowest at the trench. 

One drawback to the original wave height subroutine is the 
computational speed. When running the program to generate an energy 
spectrum for sea and swell, the model takes approximately 45 minutes per 
frequency band. Even the monochromatic surf module is slow (4.5 - 5 

minutes) because of the number of iterations required to converge to the 
preassigned criterion. 

Least acceptable of all is the Nearshore Circulation module which 
proved to be unstable as well as requiring as much as 23 hours to 
asymptotically approach an answer. To provide fast, accurate longshore 
current information, a one dimensional model is substituted for the two 
dimensional circulation model. While some potential rip current information 
will not be provided by this simple model, the one dimensional model has the 
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WAVE HEIGHT FIELD (THETA = 100) 
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Figure 3.21. VVnvo llelglits nt X = 8 with cyclic hourulary c'oiuli tioris 
irnjx)secl (solid) compnrecl with the wnve height fiehi using the 
old bouridnry conditions (dnshed) The wave pf'nks have shifted 
to the edges of the domain ns suggested by theory. 
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LONGSHORE DISTANCE (M) 

Figure 3.22. Wave Directions at X = 8 with cyclic boundary conditions 
imposed (solid) compared with original direction 
calculations (dashed). Note close agreement with 
expected behavior (dotted). 



advantage of being fast, accurate and able to provide the most impxDrtant 
input, the longshore current distribution. 

A state of the art, one-dimensional wave height model is used as 
input to the longshore current model. The new wave height and longshore 

current subroutines are described below. 

C. ONE DIMENSIONAL WAVE HEIGHT CALCULATION 

To Increase the SSSP’s computational efficiency while maintaining 
reasonable accuracy, a one dimensional random wave height model was 
substituted for the original two dimensional monochromatic model. This 

model characterizes the transformation of the wave height probability density 
function (pdf) from offshore to the shoreline (Thornton and Guza, 1983 ). 
The model assumes that the waves are narrow banded in frequency and 
direction such that the waves are specified by a single mean frequency f and 
mean direction 9 , that the bottom contours are straight and parallel, and 
that the wave conditions are stationary. This probabilistic approach to the 
wave breaking problem was originally suggested by Collins (1970) and Battjes 
(1972). The form that is used here is based on integrating the energy flux 
balance equation. The calculated wave heights are a function of the shoaling 
conditions of the integral path from offshore to the shoreline (Battjes and 
Janssen, 1978 ). 

For straight and parallel contours, the energy flux is balanced by wave 
breaking and the frictional dissipation: 



9 (E Cg cos 6) 



D 



<C^> 



(18) 
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where : 



<Gj^> = dissipation due to wave breaking 

<G^> = frictional dissipation 

In the domain of interest, it has been shown that the frictional dissipation is 
less than 3% of the dissipation due to breaking waves (Thornton and Guza, 
1983) and will be neglected. 

The waves in the new model are described by the Rayleigh wave height 
distribution: 



p(H) 



2H rms 

H ^ 
rms 



(19) 



where : 



H = root mean sauare wave height 

rms ■* ^ 

where H is the root mean square wave height. This distribution takes 
into account the fact that there is a distribution of wave heights throughout 
the entire field for a given frequency. Previous random wave models 
truncated the Rayleigh distribution in various ways. Recent observations 
have found that even after the waves break, the Rayleigh distribution- is still 
valid (Thornton and Guza, 1983). Therefore, in the improved model, the 
waves are everywhere Rayleigh distributed. The probability of a wave 
breaking described by a simple weighting, of the Rayleigh 

distribution : 

Pj^(H) = W(H) p(H) 
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The weighting function, W(H), which makes larger waves at a given depth 
more likely to break, is given by (Thornton and Guza , 1983): 



H 4 

W(H) = 



( 20 ) 



where : 



Y = adjustable parameter based on 
field data 

The model assumes that breaking waves are similar to periodic bores 
(LeMehaute', 1962) as depicted in Figure 3.23. The bore dissipation per 
unit area is calculated: 



^bore 



1 

J P9 



(h2-hi) 

^ 1^2 



3 

-Q 




( 21 ) 



where : 

Q = volume discharge across bore 

B = breaker coefficient (described 
from data) 



The volume discharge parameter is described for a linear periodic bore 
(Hwang and Divoky , 1970) by: 

Q = Ch/L (22) 

where C is the wave celerity. Using linear wave theory, the wave energy 
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-J 
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Figure 3,23. Periodic bore used to describe spillliig breakers (after 
Thornton and Guza, 1983). B is the breaker coefficient 
wliich describes how much of the wave is breaking at any 
given time. 



and group velocities are given by: 



E = p g / p (H) dH 

0 



1 

8 pg 




( 23 ) 



and , 



C 

8x 




2 k h 

Sinn 2):n 



] 



COS 9 



(24) 



After integration and the substitution of the probability density function , 
Equation 21 leads to: 









^ ^ 4,5 rms 

y h 



(25) 



where : 



f = mean frequency 

This equation for the bore dissipation can be tuned to the data by varying 
the breaker coefficient B , and the adjustable parameter , Substituting the 
bore dissipation function (25), the energy flux equation (23) is given by: 



3 x 4 » 9 Cg cos 6) 



3 £B^ 

16 P 9 4 5 ^rms 

> h 



(26) 



For shallow water, the group velocity Cg is approximately equal to the wave 

1/2 

celerity C which equals (gh) . Equation (26) was chosen to be used in 

the wave height module because, although a more complicated expression 
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gave a slightly better fit to the data, Equation (26) provides an analytical 
solution which allows a detailed study of the model's behavior. The 
analytical solution is obtained for normally incident waves impinging on a 
plane sloping beach and is used to compare the accuracy of the numerical 
model. The analytical solution is solved by integrating (26) for the wave 

height (Thornton and Guza , 1983); 



H. 



rms 



^ 1 / 5 , 9/10 
a n 



[1 - h 



23/4 



23/4 



v5/2 



) ] 



- 1/5 



(27) 



where : 



a 

tan 6 




tan 6 

f 

slope of beach 
o o 



The numerical model uses a modified Euler integration scheme to 
integrate the energy flux equation. A weighted predictor/ corrector loop is 
used to increase accuracy and stability as the water depth becomes shallow. 
The basic calculation is a forward stepping scheme with a variable grid size: 



E Cg 



2 











> 



+ <£. > ) AX 

2 1 



(28) 



Well outside the surf zone the wave height changes are very slight with 
distance as the bore dissipation is essentially zero, and are relatively 
insensitive to the grid size. However, once the waves start to break, the 
bore dissipation increases rapidly, modifying the wave heights. To increase 
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the resolution , the mesh size decreases by a factor of four when the bore 
dissipation reaches a value of 0.01. This flag value was arrived at by a 
series of model tests under different initial conditions. 

The model was tested against the analytical solution for a variety of 

initial values. The initial conditions insured that shallow water criteria were 

1/2 

met so that the shallow water approximation of Cg = (gh) , rather than 
the explicit Equation (24), could be used. The first test used a shallow 
steepness offshore wave (H /L = O.Ol) having a height of 1.0 meter 

normally impinging on a plane beach with a slope of 0.04. The period of the 

wave was 10 seconds. The numerical and the analytical solutions are 
compared in Figure 3.24. Agreement is very close through most of the 
offshore grid points. As the water depth gets very shallow, the step size is 
again reduced to insure the solution converges for the last few grid points. 

For a steep wave case (H^/L^ = 0.02), an incident wave of 2.5 meters 
was used as the offshore input (Fig. 3.24). The wave period, beach slope 
and normal incidence all remained the same. The numeric and analytic 

solutions are again in good agreement, although the numeric solution slightly 
underpredicts the wave height as it did in the case where equals 1.0 m. 
It is of interest to note that the high wave steepness case shows no growth 
in wave height due to shoaling prior to breaking. The wave heights decay 
monotonically toward the shoreline. 

Having satisfactorily compared the model results with the analytic 

solution, further improvements to the model's computational efficiency were 
made. Since the model is one dimensional, the depth variations in the 
longshore direction are averaged out. Sand bars and other longshore features 
are taken into account but offshore trenches and channels are smoothed out 
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by this process. The new height module, with this simplification, runs six 
times faster than the original module. 

The new wave height subroutine has many advantages. It uses a simple 
but robust probabilistic approach to the wave height problem. It allows 
variable bathymetry in the offshore direction and obliquely impinging wave 
fronts. It is fast, accurate, and the governing equations can easily be fine 
tuned by*' varying the parameters. 

D. LONGSHORE CURRENT CALCULATION 

The original longshore current model (Chen and Wang, 1983 ) used a two 
dimensional momentum flux balance which was numerically intensive.* The 
model converged very slowly (over 23 hours of CPU time in some cases), and 
the calculations exhibited unstable behavior (Devendorf, 1985 ). To provide 
the necessary longshore current information to an operational planner, a 
simple one dimensional model was substituted for the original two dimensional 
subroutine. The model was designed to work with the one dimensional wave 
height model discussed in the previous section. The random wave probability 
model is used as input for the longshore current calculations. A narrow 
banded wave distribution is assumed so that the wave frequencies can be 
approximated by a mean frequency f and direction 6 . The longshore 
momentum equation is a balance between the mean wave momentum flux 
divergence and the longshore component of the bed shear stress (Thornton 
and Guza, 1985): 



3S 



vx 



8X 




(29) 



71 



where : 



S 



yx 



time and depth averaged covariance 
between unsteady velocity components 



_b 

'y 



bed shear stress in longshore direction 



is made up of two terms which are assumed to be statistically 

independent. The first term is the radiation stress term (S ) which is the 

yx 

wave induced momentum flux. The second term (S is the depth 

yx 

integrated turbulent Reynold's stress, often parameterized in terms of an 
eddy viscosity term. 

Monochromatic current models require the eddy viscosity term to smooth 

out the discontinuity in the current field at the breaker line. Since no 

waves break outside the surf zone in a monochromatic model, there is no 

current generated. However, there is a jump discontinuity at the breaker 

line which requires smoothing by the eddy viscosity term. In the random 

wave model, there is no clearly defined breaker line. Some waves break 

offshore and, as the beach is approached, more waves break until, in the 

inner surf zone, almost all waves are breaking (Thornton and Guza, 1985). 

The smooth transition of energy dissipation reduces the requirement for the 

eddy viscosity term, which simplifies the velocity calculation. 

The radiation shear stress S is described by (Longuet-Higgins and 

y X 

Stewart , 1964 ) : 



S 



yx 



E C 



o 

o 



cos 6 



sin 9 

C 



(30) 
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Snell's law requires that: 



sin 0 



sin 6 



= constant = 



(31) 



so that when (30) is differentiated the result is: 



dX^ yx 



sin 6 






s in 9, 

D 



(32) 



This implies that the shear stress is balanced by the bore dissipation for 
breaking waves on a straight and parallel beach. 

The bed shear stress is given as: 



Ty = p ]U| V 



(33) 



where : 

|U| = instantaneous total bottom velocity 

V = longshore current velocity 

The bed stress is linearized by assuming small angles of wave incidence and 
weak longshore currents (Thornton, 1970; Longuet-Higgins , 1970): 



Ty = P |U| V 



(34) 



Using shallow water wave theory, the bottom velocity is computed by: 



1/2 



5 1 = |(2) I 



(35) 
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Making the appropriate substitutions into the longshore momentum relation and 



solving for V : 



V 



3 

4 



7 ^ 1/2 
B f g 



CfY 



sin 6 

C 

o 



H 



rms 

9/2 



(36) 



The numerical scheme is compared to the analytical solution, which is 

obtained ‘by simply substituting the expression for (Equation 27) into 

the velocity calculation (Equation 36). The numerical and the analytical 

solutions for the shallow and steep wave cases studied above are compared in 

Figures 3.25 and 3.26. The numerical solution of V slightly underpredicts 

the current velocity compared with the analytic solution. This is due in 

part to the introduction of a small error in H which was discussed in the 

previous section. The value of is raised to the sixth power in Equation 

(36) which means that V may be sensitive to small errors in H 

^ rms 

However, the magnitude of the error for V is on the order of 0.05 and 0.07 
m/sec which may be neglected. 

The initial runs of the updated wave height and longshore current 
subroutines are presented and the reworked SURFCON module is discussed in 
the next section. 

E. IMPROVEMENTS TO THE SURFCON MODULE 

During the reprogramming effort, a major consideration was maintaining 
the integrity of the original model, while incorporating modular 

improvements. The new height, direction and current subroutines were coded 
under new names, leaving the old subroutines intact. To maintain 

consistency, the control module, SURFCON, was renamed OLDSURFCON and 
a new module, NEWSURFCON, was written to make use of the new 
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Figure 3.25. Numeric (solid) versus nnalyUcal (dashed) loiishore 
current solutions for 11 =1.0 meter. Numerical values 

are sliglitly lower because the wave heiglil values are 
raised to the seventh power in tlie velocity calculation. 



WAVE CURRENT FIELD (THETA = 170) 

NUMERICAL VS ANALYTICAL SOLUTIONS 
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OFFSHORE DISTANCE (M) 

Figure 3.26. Numeric (solid) versus analytical (dashed) lonshore 
current solutions for II =2.5 meter. 



subroutines . 



The advantage of this approach is that the model can be 



changed from the new to the old version simply by changing the main 

program "CALL NEWSURFCON" to "CALL OLDSURFCON". 

The ISTEWSURFC ON module is basically the same as the old module with 
some exceptions. The wavenumber k, celerity C and group speed Cg are 
initially calculated. Then the new wave direction and height subroutines are 
called, which return direction, height and longshore current matrices. 

The first major change incorporated in this new module is the way the 
breaker line is defined. In the old module, a quadratic equation was solved 

for the location of the first breaker line. Since the new breaker heights are 

calculated using a weighted probability density function, P^(H) = W(H) p(H), 
there is, by definition, no specific spot where all waves are said to be 

breaking. The weighting function, W(H), as defined in Equation (37), gives 
a measure of the percent of breaking waves (Thornton and Guza , 1983 ): 



W(H) 



H 4 
r rms-| 

^ hy ^ 



= / 



0 



p^(H)d{H) 



(37) 



To specify a "breaker line" , the model first differentiates between steep 
waves, which show no wave growth due to shoaling towards the shore, and 
shallow waves, which increase in wave height until they break. The critical 
wave steepness, determined from Equation (27), is H^/L^ = 0.02. A wave 
with a steepness greater than this criteria is considered steep while a wave 
with steepness less than 0,02 is considered shallow. For steep waves, a 
decision must be made as to what percentage of waves having broken 
constitutes the surf line. After some rather ad hoc experimentation, a value 
of 0.33 was chosen as a reasonable, w^hich in effect defines a "significant 
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breaker height". That is, once 33 percent of the waves are breaking, the 
model chooses that point along the grid as the first breaker line . This test 
is actually conducted in the wave height subroutine. As the model moves 
down the column, the parameter W(H) is tested at each grid point to see if 
it greater than or equal to 0.33. Once this value is reached, the breaker 
height, breaker depth, and distance from shore are returned to the 
NEWSURFCON module. For shallow waves, the model selects the maximum 
wave height as the breaker line. These two methods give good estimates of 
the observed breaker line, as discussed in Chapter 5. 

Once the initial breaker height is defined in this manner, the number of 
additional lines are computed by dividing the initial breaker distance by the 
wave length corresponding to the average surf zone depth. If there are one 
or more wave lengths between the initial breaker line and the beach, the 
number of breaker lines are incremented accordingly. No additional 
calculations are performed on these interior breaker lines, if present. The 
breaker type is calculated based on the value of the surf parameter. Eta, 
which is a function of the beach slope and the deep water wave steepness 
(Battjes, 1974): 



= 



tan B 



tan 6 



1/2 



H 



(2tt (■ 



o,,l/2 



gT 



(38) 



If Eta is less than or equal to 0.4, the breaker type is defined as spilling. 
Eta values between 0.4 and 2.0 are defined as plunging and values greater 
than 2.0 are called collapsing or surging waves. 
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The next calculation performed by NEWSURFCON is the "Effective surf" 
computation. Effective surf is a measure of surf intensity and is used by 
operational planners to provide a criterion for the feasibility of conducting 
amphibious landing operations using specific kinds of equipment. Effective 
surf, expressed as a wave height in feet, is a modification of significant 
breaker height, taking into account the total beach and surf conditions. It is 
important ^o note that effective surf is simply a planning parameter, and not 
a direct correspondence to the significant breaker height. Most commonly, 
^eff larger than H^, but in some cases, especially when the is small 
and the breaking waves are all spilling, can be smaller than 

The calculation of effective surf uses a simple "look-up table" method 
which is outlined in the Joint Surf Manual ( COMNA VSURFPA C/ 
COMNAVSURFLANT Inst. 3840.1, 1976). The calculations of effective surf 
are retained from the original model with only minor modifications. 
Originally, the effective surf was retained as a wave height matrix. 
However, is simply a single number for a specific beach, based on 

significant wave height and direction, period, wind speed and direction, and 
percentage of spilling waves. The new module calculates a single number for 
^eff outputs all the parameters as specified by the Joint Surf 

Manual. Other information of use to the planner is printed at the same 
time. The new calculation also allows for a wind input (speed and 
direction) which was not available previously. 

The final operation in NEWSURFCON is the limitation of the wave heights 
to the maximum stable wave height for a given depth and period. This 
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maximum wave height is given as (Wang and Chen, 1983): 



= 0. 067 tanh ()ch) {^) (39) 

Each generated wave height is compared against this maximum value and is 

reset to H if it exceeds the limiting wave height. The wave height, 
max 

direction,^ current and effective surf matrices are then sent to OUTSURF for 
formatted output to either the screen or to a designated printer. 

F, DEVELOPMENT OF AN OPERATIONAL USERS GUIDE 

A simple users guide was written for the new version of the SSSP (Gill, 
1985). Several guidelines were adhered to during its development. The 
manual assumes some familiarity with the HP-9845/275 and a working 
knowledge of oceanographic terms. The manual is tutorial in nature. It is 
divided into two parts, for the two main modules OPEN and SURFCON. The 
manual explains each input request and appropriate examples are included 
where necessary. An effort was made to keep the manual brief and user 
friendly. The model, being extensively menu driven, is straight-forward 
enough to require little documentation, after it has been used once or twice. 
The BASIC code is extensively documented throughout. Where appropriate, 
the old program lines are left in place (commented to prevent execution), 
with the replacement lines documented to point out the reason for the 
specific change. 
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lY. DISCUSSION 



A. SENSITIVITY TEST OF BREAKER PARAMETERS 

As discussed in Chapter III, the computation of using the Rayleigh 

probability density function model, depends on two breaker parameters, 
Gamma and B. Gamma is an adjustable coefficient of 0(l) which appears to 
be strongly related to the beach slope (Sallenger and Holman, 1985 ). B is 
the breaker coefficient which represents the percentage of foam on the face 
of the breaking wave (i.e., the intensity of the breaking wave). This 
parameter is expected to be less than or equal to 1.0, since a B value of 
1.0 corresponds to a fully developed bore (a wave that has fully broken from 
crest to trough). Values of B greater than 1.0 imply that the bore 

dissipation function underestimates the wave dissipation due to breaking 
(Thornton and Guza , 1983). 

To verify the validity of the model assumptions, 24 wave data sets were 

studied to gain an understanding of the interaction of the B and Gamma 

values as used in the model. These data sets are summarized in Table 4.1. 

The data were digitized and appropriate Gamma values chosen, based on the 

slope of H /h inside the surf zone. The Gamma values ranged from 0.39 
rms 

to 0.60 and, for each data set, the specific Gamma value was held constant 
while the B values were varied. A curve fitting routine was used to 
optimize the B parameter for a given Gamma. As B was varied, the new 
model curve was compared with the real data points, and a sum of the 
square of the errors was generated. The B values were iterated until the 
error sum was minimized. Plots of the model output versus the data points 
were generated for each case. 
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A total of 24 data sets were tested against the model (See Table 4.1). 
In every case, the model accurately predicts the wave heights described by 
the data points. The first ten data sets are from a laboratory study of 

energy saturation of irregular waves during shoaling (Vincent, in press). 
These wave tank tests were conducted with a plane bottom slope of 1:30, 
with various frequency and incident wave heights. The experiments can be 
divided up into two groups based on wave steepness. The low steepness 

cases (H^/L^ < 0.02) are represented by Case 1638 (Fig. 4.1). The 
Gamma value of 0.60 was derived from the slope of H /h as described 
above. Using the curve fitting routine, an optimum B value of 0.82 was 
chosen. Note the good fit of the model curve (solid) compared with the 

data points (circles). The high steepness cases are represented by Case 

1148 (Fig. 4.2). Gamma was set equal to 0.60 and a B value of 1.2 was 
obtained. As before, there is good agreement between the model curve and 
the data points. 

To determine the optimized modeUs behavior when used with more 
difficult bathymetry, data from several wave tank experiments with offshore 
bars were tested (Battjes and Stive, in press). Case 5 (Fig. 4.3) has a 
simple offshore bar with a mean bottom slope of 1:25. Gamma was estimated 
to be .45 and a B value of 1.2 was obtained. Although the presence of the 
offshore bar (Fig. 4.4) radically changes the wave height data profile, the 
optimized model curve again faithfully predicts the wave heights. Wave tank 
test case 15 (Fig. 4.5) has a very complex bathymetry (Fig. 4.6). The 
model divides the on /offshore distance into 10 grid points with the 
intermediate depths calculated by straight linear interpolation. Therefore, 
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TAHLE 4.1. SURF PARAMETERS 



O' mOCCCO 

Ci3 OOOOOXOOr^r^ mcntn OL^r^^OO 

<r<J*-^<T<TXr>-r>.00 OOOOO OOOOO 

!j- ... 

0000030000 ooo ooooo ooooo 





o> 


o 


CO 


<T 




<r 


X 


CNI 


CN 


o 


cc 


X 


m 


<? 


X 


in 


lo 


X 


m 




o 


o 


o 


cn 




o 


ro 


<T 


04 


<3 




X 


o 


m 




3 


m 


X 


— 


< 


o 


m 


o 


o 


r- 


o 


o 


<r 


o 


in 


ro 


ro 


o 


04 


04 


<r 




o 




E- 


m 


m 


04 


rsi 


r>} 




— 




*— 


— 


O" 


— 


04 


o 




CO 


X 


X 


X 




o 


o 


o 


o 


O 


o 




3 


o 


o 


o 


o 


o 


— 


— * 


O 




3 




< 


o 


o 




o 


o 




o 






o 


in 


r>. 


in 


o 


X 


LO 


o^ 


'C 


o 


SI 


o 


o 


o 


o 


o 


o 


o 


O 


o 


o 


'C 


<T 


m 


'C 


<T 


'C 


'C 


X 


<0 










































<: 


o 


o 


o 


o 


o 


o 




o 


w' 


o 


3 


O 


O 


o 


3 




o 


o 


o 








m 


m 


X 


X 


X 


X 


X 


X 


X 


X 


3 


3 


3 


O 


o o 


3 


o 


o 




O 


3 


3 


3 




m 


X 


X 


X 


X 


X 


X 


'X 


X 


X 




LO 


OC 


o 




CO 


m 


X 


o 


3 


o 


*<r 


04 


3 




m 


X 


X 


X 


X 


X 


X 


X 


X 


X 


■*c 


04 


•<r 


m 




X 


X 


X 


OsJ 


OJ 


OvI 


Os} 


04 


04 


o 


o 


o 


o 


o 


o 


3 






3 






3 


o 


O 


o 


o 


3 




o 


3 


o 


3 


3 


3 


X 


3 


o 


o 




o 


3 


o 


3 


3 








3 


O 


o 


o 


o 




o 


3 


3 


o 


3 


3 





m 






04 


X 


04 


o- 




04 


3 


3 


LO 


'C 


o 


ooo 


3 


O 


3 


3 


3 


o o 




in 




CO 


X 


X 


O 


o 


O 






X 


o 


X 


X 


O — * 




in 




X 


LO 










o 


3 


3 


3 


o 


3 


O 


O 




o 


o 


o 


3 















o 



E 

u 



o 



u 

o 



05 





04 


X 


o 




04 




<r 


04 


X 


LO 


3 


04 


CO 


LO 


0 


0 


X 


CO 


— H 


X 


04 


LO 


X 


X 


o o 


04 


m 


X 




X 


04 


o 


'cr 


or 


O 


o 










0 


— • 


X 




— 




o 


_ 


— 








'•T 


LO 


<• 




— 


X 


04 


o 






3 


0 


0 


3 


0 




3 


3 


o 




o 


3 


o 




o 


w 


o 


o 


o 


3 




o 


.3 






3 


0 


0 




3 


- 


3 


o 


o 




3 


o 




o 


3 




o 


3 


O 


o 


o 


0 


0 


0 


3 


0 


3 




3 


3 


3 



0 


0 


0 


3 


0 


X 




0 


0 




'4T 


LO 


3 


0 


m 


LO 




0 


0 


LO 


r>^ 


<0 




0 












CO 


X 


X 


LO 


in 


04 


3 


X 




3 


3 


X 






04 


04 


3 


3 





0 


0 


0 


0 


0 


0 


0 


0 


'C 


<r 


X 


X 




04 


lo 


LO 


lo 


0 


04 


c^ 


X 


X 


-v^ 


04 


0 


0 


0 


3 


0 




04 


Osl 


X 


X 


X 


X 


m 


X 


X 


X 0 


04 


X 


LO 


X 


r>- 


0 * 


'"'O 




























0 


— 


— « 


lO 


0 


CO 


— 


0 


04 


04 


c^ 




























X 


X 


X 


04 


— 


X 


m 


04 






X 



X 


X 


0 - 


<r 


LO 


0 — 




lD 


0 




0 


04 


0 


0 


0 


0 


3 


0 


3 


3 


3 




3 




CO 


in 


04 




X — 






0 


04 


X 


X 






04 




0 


LO 


0 


X 


X 


c^. 


3 








— 






— N 




—I 








X 


-3“ 


LO 




04 


X 


LO 


X 


X 


"C 


lO 



0 


0 


X 


<3“ 


X 


X 0 


04 


X 0 


X 


0 




04 


X 


<r 


lo 






3 


04 




X 


■*> 


0 


04 


X 


0 


04 


LO »— 


X 


'<r 0 




— 




















— < 




04 


m 


m 


0 


0 


vC 






— 04 


V 


0 


CJ 












> 


> 


> 


> 


> 


> 




















cn 


X 


0 


c 


0 


c; 


c; 


c 


c 


c 




c 


0 


















c 




rs 


r. 


11. 


II, 




Ll, 


2 : 


zz 


zz 


zz 


zz 


2 : 



83 



(>()') U.UOl 




(WD) IHOiaH 3AVAV SWH 



84 



200 -too 600 800 1000 1200 IJOO 1600 1000 

OFFSHOliE DISTANCE (CM) 

Figure Optimized model output for Case 1638 (Vincent, 1985). 
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Figure 4,2, Oj)Llinize(l model ouLpiit for Case 1148 (Vincent, 1985). 
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OFFSHORE DISTANCE (CM) 

Figure 4.3. Optimized model output for Case 5 (Battjes and Stive, 

1985). 



Q. 

CU 

u 

oo 

o 

a. 

o 

H 

E 

O 



o 

CO 



IT) 

W 

cc 

o 




ZJ 

u 

c 



in 

00 



o 

> 



CO 



o 



XJ 

cc 

CO 



o 

E 

5^ 






CD 

> 

cc 

:i 

in 

CD 

tn 

CD 

o 



<■ 

CD 

u 

3 

00 



O O O Q 
(01) i[:jcf9g 



87 




(WD) 1H0I3H 3AVM SWH 



88 



400 000 1200 1600 2000 2100 2800 3200 3600 4000 

OrrSIlORE DISTANCE (CM) 

Figure 4.5. Optimized model output for Case 15 (Battjes and Stive, 

■ 1985). 
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Figure 4.6. Case 15 wave tank hatliyinetry . (Rattjes and Stive, 1985). 



complex bathymetry, such as that in Case 15, was smoothed out to a 
certain degree. Even with that limitation in mind, the optimized model 
curve appears to accurately estimate the recorded wave heights. The model 
results will* be tested against field data in the next chapter. 

It is evident that the random wave model, with the appropriate choice 
of the surf parameters, Gamma and B, will accurately forecast wave heights 
throughout- the domain. A sensitivity test for the choice of parameters was 
conducted on three representative cases 1638, 1148, and 5. Case 5 (Battjes 
and Stive, in press) was chosen for discussion here because of its more 
interesting bottom profile. It was felt that this barred profile was 
representative of many real beaches and would provide a more rigorous test 
than a simple plane bottom. Using a Gamma value of 0.45, the optimum B 
value of 1.3 was changed to 1.7, 1.5, 1.1 and 0.9 (Fig. 4.7). B values 
less than optimal tend to over-estimate the wave heights and vice-versa. It 
is of interest to note that changes in B of 22-33 percent only results in an 
increase in the model error of about 10-12%. This was confirmed in a study 
of Torrey Pines, California wave data by Thornton and Guza (1983). They 
found that a variation of +/- 25% about the optimum B value resulted in an 
increased model error of less than 10% (Fig. 4.8). 

It was found that the model demonstrated a similar sensitivity to changes 
in the Gamma value. A change of about +/- 20% to the optimum Gamma 
value of 0.45 resulted in an increased model error of 8-10% (Fig. 4.9). 
Sensitivity tests were conducted on high and low steepness plane beach cases 
(1638 and 1148) with similar results. In the next section, techniques and 
rationale for choosing B and Gamma values, based on initial conditions, are 
discussed. It is important to bear in mind the above sensitivity analysis. 
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Figure 4.7. SenslLivity test; G = 0.45; 

about the optinuim B value of 1. 
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Figure 4.8. Error percentage only increases about 10% as B is varied 
by 25-30% (Thornton and Guza, 1983). 
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Figure 4.9. SensitiviLy test; B = 1.3; Cania vanes from 0.35 to 
0.55 about the optimum Oama value of 0.45. The model 
ap|)ears to be slightly more sensitive to flucuations in 



The model appears to be fairly robust with respect to small changes in the 
surf parameters. Therefore, slight errors in the choice of these values 
should not degrade the final performance of the model to a great extent. 

B. CHOOSING THE BREAKER PARAMETERS 

The two breaker parameters, B and Gamma, are multiplicative constants 
in the wave height calculation (see Equation 27). They could 

therefore^ be combined as a single bulk modulus coefficient, or kept as two 
separate parameters. It was felt that keeping the two terms separate would 
better serve the requirements of the model. 

Some work has already been done on the problem of choosing the breaker 
parameters. A recent study by Sallenger and Holman (1985) has determined 
a functional relationship between Gamma and the average beach slope based 
solely on field observations inside the surf zone . This relationship is given 
by: 



Y = 3.2 tan 6 + 0.3 



(40) 



Using the 24 wave data sets, the chosen Gamma values are plotted against 
beach slope (Fig. 4.10). The solid line is a straight linear regression of the 
data points while the dashed line is the Holman/ Sallenger function. There 
seems to be a general agreement in the slope of the line although the 
intercept is slightly different. Gffmma is also plotted against wave steepness 
(Fig. 4.11) but the plot is less linearly correlated. As an initial trial, the 
Holman /Sallenger functional relationship will be adopted as the models 
convention for choosing the Gamma value. The average slope is defined here 
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GAMMA VERSUS SLOPE 
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BEACH SLOPE (Lan BgIq) 

Figure 4.10. Gamma versus bcncli slope. Solid lino is a linear 
regression fit to the data. Note tliat the single point atjovo the 
solid lino Is heavily weighted since It represents 10 cases. Dashed 
line is the llne.ar relation suggested Sallenger and Holman (108S). 
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STtEPNESS (llo/Lo) 



as the depth at the outermost grid point divided by the offshore grid 
distance . 

Choosing a proper predictor for the B value proved to be more difficult. 
Since B represents a measure of breaker intensity, it is felt that it would be 
related to incident wave steepness which has previously been shown to be 
related to breaker type (Battjes, 1972). The optimum B values for the 24 
wave sefs are plotted against steepness in Figure 4.12 and against beach 
slope in Figure 4.13. The wave steepness in this case was defined as the 
incident wave height divided by the deep water wave length (H^/L^). 
Although there is a considerable scatter, B appears correlated with wave 
steepness and less with beach slope. 

Several other parameters were considered as predictors for B. The B 

and Gamma values were compared with the surf parameter, Eta = tan B / 

1/2 3 4 

(H^/L^) . Additionally, the Bulk Modulus (B /Gamma ) was plotted 

against beach slope and wave steepness. Since, from Figures 4.11 and 4.12, 

B appears to be related to both beach slope and wave steepness, it was 

hoped that Eta, which is a measure of the expected breaker intensity, would 

offer a predictor for B. The plot of B versus Eta (Fig. 4.14) shows an 

interesting pattern. B appears to increase almost linearly as Eta increases 

to about 0.6 and then decreases in an exponential fashion to an asymptotic B 

value of 0.8. This behavior appears reasonable. As the surf parameter 

increases, indicating a transition from spilling to plunging waves, the B value 

(amount of wave face that is breaking), should increase. 

The maximum value that B should reach, according to theory, is 1.0. 
Studies have shown, however, that in many cases, the optimum B value is 
on the order of 1.4 (Thornton and Guza, 1983). The maximum B value 
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BEACH SLOPE (tan Bota) 
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Figure 4.13. B ai)()ear.s Lo Be a function of incideuL wave steepness. 
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SURF PARAMETER (ETA) 

Figure 4.14. B versus the surf parameter Eta. B is strongly related to 
Eta which is to be expected. Both B and Eta are measures of 
tlie intensity of the wave breaking. Dashed line depicts 
piecewise function used by the model to define B. 
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SURE PARAMETER (EIA) 

Figure 4.15. Gamma versus Fta . 



occurs when Eta is about 0.6, indicating that the breakers are mixed spilling 
and plunging. As Eta continues to increase, the B values fall more into line 
with theory and approach the expected value of 1.0, as the waves become 
fully plunging. 

To simplify and expedite the model development, the following piecewise 
continuous function is used to determine the B term (Fig. 4.14): 

B = 2.167 Eta + 0.4 Eta <= 0.6 (4l) 

B = - 1.70 Eta + 2.7 0.6 < Eta <= 1.0 

B = 1.0 Eta > 1.0 

This choice is somewhat ad hoc and is based primarily on the data sets 
discussed above. Future refinements in the techniques for predicting the 
breaker parameters may substantially improve the use of this model 
operationally. Using equations (40) and (41) as predictors for Gamma and 

B, the SURFCON module is essentially complete. The next chapter will 
discuss testing the model against two well-documented natural beach data 
sets, to get a feeling for how the improved model behaves against actual 
data. 

C. MOTIVATION FOR THE RANDOM WAVE MODEL SUBSTITUTION 

As discussed in Chapters 3 and 4, the primary motivation for substituting 
the one dimensional random wave model for the two dimensional model was to 
increase computational efficiency and accuracy. The improved model 
calculates longshore current information quickly and accurately. Significant 
wave height predictions compare well with actual data sets, as will be 
discussed in Chapter 5. There is , however, a potential for loss of detailed 
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information by using the one dimensional model. Nearshore circulation aind 
longshore variation in wave heights are smoothed out by the one dimensional 
approach , 

The surf model (both old and new versions) uses a single wave height 
input to define the offshore boundary condition. This implies that variations 

in the interior wave heights and directions are due only to variations in 

bathymetry. Because the model uses a very small grid (only 13 longshore 
grid points), very accurate bathymetry must be entered. It is questionable 
whether the operational planner would have access to bathymetric data of 
the required accuracy to define a two dimensional bathymetry. 

Many beaches which would potentially be used for amphibious operations 
are shallow, almost planar in nature. Other common landing beaches are 
barred, but rather uniform in the longshore direction. Therefore, for many 
operationally relevant beaches, the assumption of straight and parallel 
contours, and the use of a one dimensional model, appears justified. The 
"mean beach" surf information provided by the model will adequately 
represent the conditions over a standard 500 yard landing beach. If a wider 
beach needs to be characterized, the improved model would be run several 

times, entering the required information for each beach segment. 
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V. TEST AND EVALUATION OF THE IMPROVED MODEL 



A. WAVE HEIGHT AND BREAKER LOCATION 

The improved version of the SSSP was tested against one wave tank 
experiment and two sets of beach surf data to see if the model would choose 
reasonable^ Gamma and B parameter values, based on the predictors discussed 
in Chapter 4, The predicted wave heights and surf characteristics are 
compared with the observed values. 

Case 5 (Battjes and Stive, 1985), discussed in the preceding chapter, 
was chosen as the wave tank experiment with which to test the finished 
model. This case was chosen because of its barred beach bathymetry. 
The model selected values of 0.41 and 0.963 for Gamma and B. The 
optimum values, from Table 4.1, are 0.45 and 1.30. The model wave height 
fit to the observed data points is still quite reasonable (Fig. 5.1). The 
solid line represents the model output and the circles are the observed data. 

Two Southern California beach data sets were used for further tests. 
Data collected at Santa Barbara's Leadbetter Beach during February, 1980 
and Torrey Pines during November, 1978 (Thornton and Guza, 1983), were 
used as the primary tests of the performance of SSSP. 

The first Santa Barbara test used data collected on 3 February, 1980. 
The average beach slope was 0.047 and the incident waves were 0.55 meters 
(rms) with a period of 14.3 seconds. The observed breaker height was 
measured at 0.70 meters, approximately 40 meters from the beach (dashed 
line). The wave angle was 7.8 degrees from normal. The SSSP predicted a 
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MODEL OUTPUT VS WAVE TANK DATA 
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Figure 5.1. Compari.son between the model output (solid) and the 
observed wave tank data (daslied). C, for all cases is 
0.009. ' 



rms breaker height of 0.68 meters (solid line), 40.0 meters from the beach 
(Fig. 5.2). The model under predicted the observed heights less than 4 
percent . 

Similar results are obtained for Leadbetter Beach data obtained on 4 and 
5 February, 1980. On 4 February, the incident wave height was 0.56 

meters with a period of 14.3 seconds. The wave angle was 9 degrees from 
normal. The observed breaker height was 0.67 meters, 50 meters from the 
beach. The model predicted 0.65 meters, 40 meters from the shore (Fig. 
5.3). 

On 5 February, the incident wave was 0.45 meters with a period of 
12.8 seconds. The wave angle was 8.4 degrees from normal. Observed 

breakers were 0.60 meters, about 36 meters from the beach. SSSP predicted 
0.55 meter breakers, 33 meters from the shore (Fig. 5.4). 

Torrey Pines beach data acquired on 4 and 10 November, 1978 is 
compared with model results in Figures 5.5 and 5.6. The beach has a shallow 
slope of 0.022 which accounted for a wide surf zone of about 160 meters. 
On 4 November, the incident waves were about 0.35 meters with a period of 
14.3 seconds and normal incidence. The observed breaker height was 0.47 
meters about 110 meters from the shore. SSSP predicted 0.48 meter 
breakers at a distance of 106 meters from the beach. The model correctly 
forecast the wave height and surf zone width (Fig. 5.5). 

A similar result is seen for the November 10 Torrey Pines Case. The 
incident wave is 0.68 meters with a period of 15.9 seconds and normal 
incidence. Observed breakers were 0.72 meters, approximately 130 meters 
from the shoreline. The model predicted breaker heights of 0.75 meters, 140 
meters from the beach (Fig. 5.6). As for the November 4 case, the model 
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Figure 5.2. Comparison between the model output (solid) and the 
observed data from Santa Barbara on 3 Fob., 1^)80 
(dashed) . 
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Figure 5,3. Comparison between the model output (solid) and the 
observed data from Santa Barbara on 4 Feb., 1980 
(dashed) . 
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observed data from Santa Barbara on 5 Feb., 1980 
(dashed) . 
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Figure 5.5. Comparison between the model output (solid) and the 
observed data from Torrey Pines on 4 Nov., 1978 
(dashed) . 
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Figure 5.6. Com|)arison between the model output (solid) and the 
observed data from Toney Pines on 10 Nov., 1978 
(dashed) . 



worked quite well for wave height prediction and the location of the breaker 
line . 

Based on the tests of the SSSP against various data sets, it appears that 
the new version is an improvement over the original model. The old version 
consistently underpredicted the wave heights and was very inaccurate at 
estimating the location of the breaker line (Devendorf, 1985). The new 
version appears to be quite good at predicting the breaker height and gives 
better estimates of the breaker line. Table 5.1 is a comparison between the 
optimum values of the two breaker parameters and the values chosen by the 
model. The technique for choosing Gamma and B is convenient, although 
improvements will need to be made in this area. 

B. LONGSHORE CURRENT DISTRIBUTION 

Longshore currents are the result of waves impinging at an oblique angle 
to the shoreline. All the wave tank tests, as well as the Torrey Pines 
cases, had normally incident waves. The Santa Barbara Leadbetter Beach 
data is the only comprehensive field data on longshore currents acquired to 
date. There exists no laboratory data for longshore currents generated by 
obliquely incident, random waves. The same three data sets (February 
3,4,5) were chosen as study cases because they had a variety of incident 
wave heights and wave angles. For the February 3 case (Fig. 5.7), the 
maximum observed current (dashed) was 0.48 m/sec at a distance of 22 
meters from the shoreline . The bed shear stress coefficient , C^, was set 
equal to 0.009 for all cases. This value was derived experimentally by 

Thornton and Guza (1985). The SSSP prediction (solid) for longshore current 
was a maximum of 0.5 m/sec, 25 meters from the beach (solid line). 
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For the February 4 case (Fig* 5.8), the maximum observed current was 
.48 m/sec, 25 meters from shore; the model prediction was .44 m/sec at 
20 meters from the beach. The Feb 5 case (Fig. 5.9) shows an observed 
current maximum of .36 m/sec at 20 meters from the beach. The model 

predicted .36 m/sec at 20 meters from the shoreline. 

The model not only predicts the maximum current velocity reasonably 

well, but ‘^the current distribution within the surf zone shows good agreement 
with the data. The model seems to underpredict the current velocity 
seaward of the maximum current, but, in general, the model does quite a 

good job of longshore current prediction. 

The other Leadbetter Beach cases showed similar results. In all cases, 
the model gives reasonable maximum current locations and speeds as well as 
current distributions within the surf zones. These findings were also 

confirmed by Thornton and Guza (1985) in a study in which the random wave 
model’s wave height and longshore current predictions are rigorously tested 
against the Leadbetter Beach data. 



114 




(D3S/WD) INaHHnO 



115 



Figure 5.9. Comparison between the model output (solid) and tlio 
observed data from Santa Barbara on 5 Feb., 1980 
( triangles) . 



MODEL OUTPUT VS OBSERVED DATA 

3 FEB 1980 B = i.OO GAMMA = .45 
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Figure 5.7. Comparison between the model output (solid) and the 
observed data from Santa Barbara on 3 Feb., 1980 
( triangles) . 
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Figure 5.8. Comparison between the model output (solid) and the 
observed data from Santa Barbara on 4 Feb., I960 
( triangles) . 



VI. SUMMARY AND CONCLUSIONS 



A. SUMMARY 

In 1981, the U.S. Navy contracted to have a Sea, Swell and Surf 
Program (SSSP) developed to provide input to the planning efforts of Naval 
Oceanographers at sea and amphibious operational planners. The model was 
developed in FORTRAN and designed to run on a main frame computer (Wang 
and Chen, 1983). The completed SSSP was then translated into BASIC for 
use on the currently available shipboard micro-computers, the HP-9845B/275 
(Devendorf, 1985). The SSSP was modular in nature, with separate 

subroutines for offshore wind wave generation , surf zone calculations and 
nearshore circulation. 

The original SSSP had several problems which prevented its use by 

operational planners. The offshore wind wave generation module overbuilt 

the wave heights for high wind speeds. Some of the surf condition 
calculations were inaccurate and the nearshore circulation subroutine had 
stability problems. These problems were addressed by this thesis in an 

attempt to improve and upgrade the SSSP. 

The open water module shows satisfactory wind wave generation for wind 
speeds of 15 knots or less. For higher wind speeds, the predicted wave 
generation occurs too quickly, both spatially and temporally. Sensitivity 
analysis indicated that the Phillips linear term was responsible for the rapid 
development, but no improvements to the model physics were made. Several 
ad hoc "fixes" were recommended. This problem will require further 

research . 



118 



Two subroutines are called to calculate the wave transformation as the 



wave proceeds from offshore to onshore. Several numeric improvements were 
made which eliminated minor instabilities in the wave direction and height 
subroutines. The original no flux boundary conditions were changed to cyclic 
boundary conditions in both subroutines. This change allows the wave 
heights and directions to more accurately reflect the conditions required by 
the model physics at the model boundaries. 

The most serious problem, the instability of the nearshore circulation 
module, is solved by substituting a one dimensional longshore current model 
for the original two dimensional scheme. Although there is a loss of two 
dimensional current resolution within the surf zone due to the assumptions of 
straight and parallel contours, stationary wave conditions, and shallow water 
initial waves, it was felt that a model that gives the important longshore 
current information quickly and accurately is an improvement over the 
original unstable, two-dimensional model. The new current model is a 
simple, longshore momentum flux balance equation. It retains the wave 
induced momentum flux and neglects the turbulent Reynold's stress. The new 
model is fast and accurate , requiring only a few calculations once the wave 
height field is generated. 

To take advantage of state of the art breaker height modeling, a random 
wave model was implemented for the calculations of the wave heights in the 
nearshore area. The random wave model assumes that wave heights, in the 
nearshore region, follow a Rayleigh distribution and that, even after the 
waves break, this distribution is retained. As a future project, the random 
wave model could be developed into a two dimensional model by applying the 
proper longshore numerical differencing. 
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The random wave model assumes that wave breaking can be described by 
a periodic bore. Two parameters are introduced by the model, which must 
be chosen based on the initial conditions. The Gamma term is part of the 
proportionality coefficient which modifies the initial Rayleigh distribution, and 
is found to be strongly correlated with the beach slope. The breaker 
parameter, B, describes the amount of the. wave face that is breaking and 
appears to be related to the surf parameter Eta. The criterion for choosing 
the two breaker parameters was developed empirically. 

The wave height and current models were next tested against wave tank 
and natural beach data sets. Although the model did not necessarily choose 
optimum values for the breaker parameters, the model is fairly robust, and 
slight deviations in the parameter values do not degrade the wave height 
calculations significantly. Good agreement was found between the observed 
wave heights and the model’s predicted values. The longshore current, based 
on the model's predicted wave height field, were also in close agreement 
with observations. 

The calculation of Effective Surf, an amphibious planning guideline, was 
improved. The algorithm was implemented in accordance with the Joint Surf 
Manual (1976) with respect to its input requirements and its formatted 
output . 

A simple users manual, tutorial in nature, was developed for this version 
of the SSSP. It is short, and assumes some knowledge of the host computer 
and basic oceanographic terms. It is written within the guidelines set forth 
by the Navy for development of computer models for the Shipboard Numerical 
Aid Program (Brown, 1984). 
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B. CONCLUSIONS 



The improved Sea, Swell and Surf Program (SSSP) can provide much 
useful information to the Naval Oceanographer at sea or amphibious planner. 
It has several problems which have not been resolved, so its use is 
recommended with several caveats. The offshore wind wave generation 

module is felt to be accurate for low wind speeds only. Useful information 
on offshore wave growth can be obtained for input wind speeds of 15 knots or 
less. The results obtained for higher wind speed inputs must be examined 
carefully before use . The requirement for providing an initial wave energy 
spectrum is probably impractical for shipboard use. The option of providing a 
single initial wave height and period will be of more value to the operational 
user. 

The new random wave breaker height model appears to work well. The 

breaker heights and locations as well as the longshore current field are 

accurately predicted by the model. Due to its one-dimensional nature and 
the assumption of straight and parallel contours, longshore wave height and 
current resolution is restricted. However, the assumption that an amphibious 
landing beach is uniform in the longshore direction is a reasonable one over a 
wide range of operational areas. The small grid size of the model supports 
the use of a one dimensional model to compute an average set of charac- 
teristics for a "mean" beach. The speed suid stability of the new model 

make it a viable tool for operational planners. Recommendations for further 
improvement of the SSSP include: 

1. Pursue the problem of the overbuilding of the offshore wind waves. 

2. Development of better predictors for the breaker parameters Gamma 
and B by including a wider data base. As was seen in the 
optimized tests, correct choice of these parameters will provide a 
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very good fit to observations, even with very complicated 
bathymetry. 

3. Explore the feasibility of using the random wave model two 
dimensionally. This would offer the advantages of this more realistic 
approach to wave height modeling, while retaining the details that 
two dimensional spatial resolution provides. 

The improved SSSP offers the Naval Oceanographer at sea and amphibious 
operational planners a new and useful tool that, when used in conjunction 
with surf observations, weather information, and other available data, can be 
a great aid in the safe planning and execution of beach landing operations. 
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